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An  Absorbing- Generating  Boundary  Condition 
for  Shallow  Water  Models 


A.R.  Van  Dongeren^  and  I.A.  Svendsen,^  Member,  ASCE 

ABSTRACT:  An  absorbing-generating  boundary  condition  based  on  the 
Method  of  Characteristics  is  derived  for  the  2D-horizontal  nonlinear  shal¬ 
low  equations.  It  allows  outgoing  waves  to  leave  the  computational  domain 
through  the  boundaries  with  a  minimum  of  reflection  while  specifying  in¬ 
coming  waves  at  the  same  boundaries.  The  boundary  condition’s  absorbing 
properties  are  tested  for  both  linear  and  nonlinear  waves  for  a  range  of  angles 
of  incidence.  Its  performance  is  compared  to  the  classical  Sommerfeld  radia¬ 
tion  condition  for  the  linear  case.  Also,  a  case  of  simultaneous  absorption  and 
generation  of  waves  at  the  same  boundary  is  analyzed. 


1  Introduction 


When  analyzing  nearshore  problems  using  numerical  models  it  is  usually  nec¬ 
essary  to  limit  the  computations  to  a  small  region  around  the  area  of  immediate 
interest.  This  implies  introducing  artificial  boundaries  for  the  computational  re¬ 
gion  that  form  the  interface  to  the  exterior  which  is  not  modelled  or  is  modelled 
in  a  simplified  way.  Thus  one  of  the  most  important  problems  in  developing  time- 
dependent  shallow  water  models  is  the  specification  of  accurate  boundary  conditions 
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along  these  artificial  boundaries  because,  after  long  enough  time,  the  performance  of 
these  conditions  will  dominate  the  model  results  in  the  entire  computational  domain. 

These  time-dependent  models  essentially  solve  approximations  to  the  equations 
of  motion  (conservation  of  mass  and  momentum)  by  integration  in  time  of  cer¬ 
tain  dependent  variables,  typically  surface  elevation  and  horizontal  velocities.  The 
boundary  conditions  are  required  to  provide  a  similar  upgrading  in  time  of  the  same 
variables  along  the  boundaries. 

In  developing  the  nearshore  circulation  model  SHORECIRC  (Van  Dongeren  et 
al,  1994,  1995),  which  is  capable  of  describing  a  number  of  phenomena  (such  as 
edge-waves,  surfbeat,  longshore  currents,  shear  waves,  etc.),  we  encountered  exactly 
this  problem.  For  our  purposes  it  was  necessary  to  develop  boundary  conditions  on 
the  artificial  boundaries  that  are  able  to  generate  a  specified  long  wave  and  simul¬ 
taneously  absorb  outgoing  waves,  i.e.  an  absorbing-generating  boundary  condition. 

Most  of  the  existing  literature  on  the  topic  of  artificial  boundaries  is  concerned 
with  absorbing  (sometimes  called  radiating,  non-reflective  or  open)  boundary  condi¬ 
tions  specifically  derived  for  the  wave  equation  or  the  shallow  water  equations.  For 
a  thorough  review  on  this  subject,  we  refer  to  Givoli  (1991). 

In  one  of  the  most  frequently  quoted  papers,  Fngquist  &  Majda  (1977)  [F&M] 
developed  a  perfectly  absorbing  boundary  condition  which  is  nonlocal  in  space  and 
time.  This  means  that  the  complete  time  history  along  the  entire  boundary  is  re¬ 
quired  in  order  to  update  the  variables  at  any  point  along  the  boundary  in  time. 
Because  this  is  very  impractical  for  any  numerical  application,  F&M  derived  local 
approximations  to  the  general  solution  of  increasing  order  of  accuracy.  The  approx¬ 
imations  are  centered  around  chosen  angles  of  incidence  On  between  the  boundary 


normal  and  the  direction  of  the  outgoing  waves.  Higdon  (1986,  1987)  derived  a 


general  form  of  this  radiation  condition  and  showed  that  it  can  be  written  as 

d  ^  Co  d 
dt  cos  On  dx 

where  n  is  the  order  of  accuracy  and  Co  is  the  linear  phase  speed.  This  expression 
gives  the  best  absorption  when  the  angle  of  the  outgoing  wave  to  the  normal  is 
For  On  =  0°  the  equation  reduces  to  E&M’s  boundary  condition,  which  only  absorbs 
normally  incident  waves  optimally.  To  the  first  order  of  the  approximation  (n  =  1), 
the  E&M  boundary  condition  further  reduces  to  the  Sommerfeld  radiation  condition 
(Sommerfeld,  1964) 


(2) 


which  essentially  states  that  the  outgoing  wave  is  propagating  in  the  positive  direc¬ 
tion  without  change  of  form.  Eq.  (1)  was  also  found  independently  by  Keys  (1985) 
and  is  an  improvement  over  the  boundary  condition  developed  by  E&M  because  the 
reflection  coefficient  can  be  greatly  reduced  if  the  angle  of  incidence  On  is  known 
in  advance.  This  might  be  the  case  for  some  types  of  problems  (e.g.,  waves  radiat¬ 
ing  from  a  source  inside  the  computational  domain),  but  not  for  the  more  general 
models  we  are  considering  here. 


Another  disadvantage  of  this  type  of  boundary  condition  is  that  the  solution  to 
the  problem  is  assumed  to  have  a  certain  form.  Broeze  &  Van  Daalen  (1992)  did  not 
make  that  assumption  and  derived  a  boundary  condition  from  the  local  energy  flux 
in  the  normal  direction  to  the  boundary  using  the  variational  principle  and  showed 
improved  accuracy  when  used  in  a  panel  method. 

Unfortunately,  the  above  mentioned  boundary  conditions  are  only  capable  of 
absorbing  waves  and  (except  for  Broeze  &  Van  Daalen)  only  address  the  linear 
problem.  Hibberd  (1977)  considered  the  more  general  problem  of  simultaneously 
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generating  and  absorbing  waves  and  derived  a  boundary  condition  for  the  nonlinear 
shallow  water  equations  (NSW)  where  the  outgoing  wave  is  calculated  using  the 
Method  of  Characteristics  while  an  incoming  Riemann  variable  associated  with  an 
incident  uniform  bore  is  specified.  On  a  similar  basis  Verboom  et  al.  (1981)  gave 
more  general  expressions  for  weakly-reflective  boundary  conditions  based  on  the 
specification  of  incoming  Riemann  variables.  Verboom  &  Slob  (1984)  derived  two 
orders  of  approximation  of  this  type  of  boundary  condition  and  calculated  reflection 
coefficients  which  were  of  the  order  of  a  few  percent.  However,  the  applications  in 
both  papers  only  deal  with  situations  where  no  incoming  wave  is  specified  and  the 
boundary  condition  reduces  to  the  case  of  absorption  only. 

In  the  case  of  simultaneous  absorption  and  reflection,  it  is  not  possible  to  specify 
the  incoming  Riemann  variable  since  it  is  a  function  of  the  unknown  surface  elevation 
and  velocity.  Instead,  Kobayashi  et  al.  (1987)  used  the  outgoing  characteristic  and 
substituted  a  linear  long  wave  relationship  between  the  velocity  and  the  surface 
elevation  to  solve  for  the  outgoing  wave.  They  only  solve  the  problem  for  one 
horizontal  dimension,  equivalent  to  the  case  of  normal  incidence  and  did  not  report 
on  the  accuracy  of  the  boundary  condition.  The  present  paper  gives  an  extension  of 
Kobayashi  et  al.  (1987)’s  boundary  condition  to  the  general  case  of  two  dimensions 
and  expands  the  condition  to  a  higher  order  of  approximation. 

The  outline  of  this  paper  is  as  follows.  In  Section  2  we  discuss  the  formulation  of 
the  problem.  In  Section  3  the  boundary  condition  is  derived  from  the  fundamental 
equations  for  two  orders  of  the  approximation.  In  Section  4  the  reflection  properties 
of  both  versions  are  investigated  for  the  case  of  absorption  only  and  compared  to 
the  classical  Sommerfeld  radiation  condition  for  the  linear  case.  In  Section  5  the 
boundary  condition  is  further  tested  for  the  case  of  simultaneous  generation  and 
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absorption  at  the  same  boundary. 


2  Formulation  of  the  problem 


The  boundary  conditions  we  specify  along  the  artificial  ocean-side  boundaries 
must  guarantee  a  unique  and  well-posed  solution  to  the  differential  equations.  As 
may  also  be  inferred  from  the  literature  review  given  above,  this  is  not  a  straight¬ 
forward  problem,  and  it  appears  that  to  some  extent  waves  and  currents  need  to  be 
addressed  separately. 

From  the  outset  one  would  expect  that  the  idea  of  emulating  the  effects  of  a  large 
ocean  in  a  computation  that  only  covers  a  small  region  of  that  ocean  imposes  some 
limitations  on  what  can  actually  be  represented  in  the  model,  and  this  is  true.  More 
importantly,  however,  it  also  requires  a  clarification  of  which  physical  mechanisms 
we  should  actually  try  to  describe  along  those  boundaries. 

Our  requirements  can  be  formulated  by  stating  that  the  boundary  conditions 
need  to  satisfy  two  criteria: 

1.  The  region  outside  the  computational  domain  can  only  influence  the  motion 
inside  through  the  incident  (long)  waves  and  through  the  currents  along  the 
boundaries.  Thus  we  must  assume  that  we  know  and  can  specify  those  currents 
and  incident  waves. 

2.  (Long)  waves  propagating  out  of  the  computational  region  must  be  allowed  to 
propagate  freely  through  the  open  ocean-side  boundaries  with  minimal  reflec- 


It  turns  out  that  whereas  outgoing  waves  can  be  separated  from  the  total  signal 
and  absorbed  at  the  artificial  boundaries,  this  is  not  the  case  for  currents.  The  dis¬ 
tribution  of  the  currents  is  essentially  an  elliptic  problem  and  therefore  the  currents 
have  to  be  specified  along  the  entire  boundary  to  uniquely  specify  the  problem.  This 
also  implies,  however,  that  in  the  general  case  extensive  information  is  needed  about 
the  currents  outside  the  domain  in  order  to  be  able  to  specify  the  currents  along  the 
computational  boundary. 

Finally  it  raises  the  question  of  how  to  distinguish  between  waves  and  time- 
varying  currents.  A  closer  inspection  of  this  problem  suggests  that  this  is  a  matter 
of  the  time  scale  of  the  variations  relative  to  the  time  it  takes  for  a  disturbance  to 
propagate  across  the  computational  domain.  For  simplicity,  we  limit  our  scope  in 
the  present  paper  to  the  case  of  incident  long  waves  without  currents. 

Thus  the  boundary  conditions  must  be  able  to  generate  a  specified  long  wave 
and  simultaneously  absorb  outgoing  waves,  in  the  presence  of  known  currents  and 
ideally  without  much  additional  computational  effort. 

3  Derivation  of  the  boundary  condition 


The  governing  equations  of  the  SHORECIRC  model  are  the  depth-integrated, 
shortwave  averaged  continuity  and  momentum  equations  (Van  Dongeren  tt  1994, 
Svendsen  &  Putrevu,  1995).  If  we  place  the  open  boundaries  carefully  we  can  achieve 
that  near  these  boundaries  the  local  forcing  is  weak.  This  means  that  the  dominating 
terms  in  the  continuity  and  momentum  equations  are  the  terms  corresponding  to 
the  nonlinear  shallow  water  (NSW)  equations  in  matrix  form: 


dt 


u 

V 

h 
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0 

9  ' 
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V 

0 

0  ' 

d 

0 

u 

0 
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+ 

0 

V 
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— 

.  ^ 

0 

u 

dx 

h 

_  0 

h 

V 

dy 

.dho 


^  oho  1  r 

9-^  +  h 
0 


(3) 


where  h  is  the  total  water  depth  h  =  K  +  C,ho  is  the  still  water  depth  and  (  is  the 
shortwave  averaged  surface  elevation,  u  and  v  are  the  depth  averaged  and  shortwave 
averaged  velocities  in  the  x  and  y  directions,  respectively.  See  Fig.  1  for  a  definition 
sketch.  It  may  be  worth  emphasizing  that  since  we  are  considering  the  general  case 
of  short  wave  motion  with  arbitrary  time  variation,  the  shortwave  averaged  (  and 
u,  V  represent  surface  elevation  and  particle  motion,  respectively,  in  the  infragravity 
wave  motion.  Usually  ^  will  also  include  a  steady  set-down  or  set-up  component,  f 
represents  all  the  local  forcing  terms  for  the  motion,  which  comprise  the  radiation 
stress  gradients,  the  current-current  and  current-wave  integrals  (originating  from 
the  non-uniform  variation  of  the  velocities  over  depth)  and  the  bottom  and  wind 
shear  stresses.  These  effects  are  all  included  in  the  original  equations. 


Following  the  procedure  outlined  in  Whitham  (1974)  and  Abbott  (1979),  the 
eigenvalues  and  eigenvectors  are  obtained  from  this  system  of  equations.  The  three 
eigenvectors  span  the  space  P: 


^  cos  cos  19  sin  i9  ^ 
sin  i9  sin  I?  —  cos  i9 

V  vf  ^  J 


(4) 


®  where  ??  is  the  angle  between  the  normal  to  the  boundary  and  the  x  axis  as  identified 

by  Verboom  et  al.  (1981)  (see  Fig.2).  Premultiplying  the  system  of  equations  (3) 
with  P~^  yields  the  governing  equations  in  characteristic  form  as  derived  by  Verboom 
®  et  al.  (1981).  (Please  note  that  the  eigenvectors  in  their  Eq.  (11.5)  contain  a  typo 

and  are  therefore  not  identical  to  the  matrix  P.) 
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It  is  convenient  to  choose  the  x-axis  normal  to  the  seaward  boundary  of  our 
rectangular  domain,  which  sets  =  0.  The  equations  in  characteristic  form  then 


simplify  to: 


dt 


-{u  -  c)— - V 


dx 


dy 


dv 


dhp 

dx 


+  Fy- 


(5) 


d/3+ 

dt 


-(«  +  c) 


d^+ 

dx 


dv  9ho 


(6) 


^7 

'm 


_d'y  _dj 


(7) 


In  (5)  the  Riemann-invariant  /?  is  defined  as 

r  =  u-2c  =  -  2\/^(^7T7)  (8) 

{ho  +  C) 

where  Qx  is  the  total  flux  in  the  x  direction  and  u  is  the  depth  averaged  velocity. 
The  Riemann-invariant  of  (6)  is  similarly  defined  &s  =  u  +  2c.  It  turns  out  that 
the  7-equation  is  the  ?/-momentum  equation  itself  which  has  the  Riemann-invariant 


7  =  n 


Qy 

{ho  +  C) 


(9) 


The  definition  sketch  in  Fig.  3  shows  that  /?■*■  propagates  along  a  characteristic 
in  the  positive  x  direction,  in  the  negative  x  direction  and  7  in  the  y  direction. 
The  forcing  terms  Fp-  and  F~^  originate  from  the  /-terms  in  (3).  These  terms 
imply  that  and  7  vary  along  their  characteristics  and  the  invariants  should 

therefore  actually  be  called  Riemann  variables. 

During  the  computation  we  will  at  time  step  n  know  the  total  value  of  (f"  and 
{Qx^  Qy)  interior  as  well  as  boundary  points.  The  incoming  wave  motion  is  spec¬ 
ified  along  the  a;  =  0  boundary  through  specification  of  {Qx,i,Qy,i)-  The  governing 
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equations  (3)  will  then  provide  the  values  of  the  total  (,  Qx,  Qy  for  interior  points  in 
the  domain,  and  the  problem  is  to  determine  the  equivalent  total  values  along  the 
boundary  at  time  step  n  +  1.  In  this  process  we  also  determine  the  parameters  of 
the  outgoing  wave. 

In  the  following,  the  absorbing-generating  boundary  condition  for  two  different 
orders  of  the  expansion  are  derived  for  a  boundary  along  the  y-axis  (x  =  0)  only. 
The  results  can  readily  be  generalized  to  arbitrary  boundary  directions  which  is 
omitted  here. 

Assuming  linear  superposition  of  the  incoming  wave  (subscripted  i  in  the  follow¬ 
ing)  and  the  outgoing  wave  (subscripted  r)  we  can  write 

Q  =  Q:c,i  +  Q.,r  C  =  (10) 

Without  further  approximation  the  incoming  Riemann-variable  (8)  can  then  be 
rewritten  as 


where 

Co  =  (12) 

At  this  point  we  introduce  the  assumption  that  the  total  volume  flux  in  the  direction 
of  the  wave  propagation  is  related  to  the  surface  elevation  by  the  equation; 

g  =  c(c  -  c')  +  Q  (13) 

where  cC  represents  the  volume  flux  in  the  oscillatory  part  of  the  motion.  Q  is  the 
net  volume  flux  which  consists  of  the  nonlinear  mass  flux,  in  the  infragravity 
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waves  and  the  ’’current”.  (  is  the  average  over  the  infragravity  wave  period  of  the 
infragravity  wave  surface  elevation  (.  It  has  been  shown  (see  e.g.  Svendsen  (1974) 
or  Svendsen  &  Justesen  (1984)  for  two  different  derivations)  that  this  relationship, 
which  is  purely  kinematical,  is  exact  for  plane  waves  of  constant  form,  no  matter 
what  height  or  nature.  Thus  the  use  of  this  relationship  here  only  implies  that 
assumption  for  the  incoming  and  outgoing  wave  motion  in  the  neighborhood  of 
the  boundary  in  question.  (Svendsen  &  Justesen  (1984)  found  that  even  for  waves 
deforming  rapidly  towards  breaking,  the  error  from  using  this  relationship  was  less 
than  5%). 

For  simplicity  we  assume  in  the  following  that  (  as  well  as  Q  are  zero.  Then, 
again  using  linear  superposition,  Eq.  (13)  for  the  x-component  of  Qi  and  Qr  can  be 
written  as 


Qx,i  =  C  0  cos  Bi  Qx,t  =  -cCt  cos  Br  (14) 

where  Bi  and  Br  are  defined  as  the  angles  between  the  normal  to  the  boundary  and 
the  incoming  and  outgoing  waves  in  the  range  [— f,  |],  respectively.  Eq.  (11)  then 
becomes 

_  Qx,i  A  Qx,i _ Qx,r  ^  Qx,r  _ 

~  Coho\^hoCCOsBi  hoCCOsBr)  Coho\  hoCCOsBi  hoCCOsBr) 


-  2 


Qx,i 

hoC  cos  Bi 


Qx,r  \  ^ 
hoC COS  Br) 


(15) 


Here  we  can  expect  that  QxlKco  <1.  If  we  expand  this  expression  to  first  order 
with  respect  to  QxjhoCo  and  solve  with  respect  to  Qx,r  we  get 

ms  H_  /  .  .\  ^  f  Qx  ''  ^ 

Qx,r 


(cos  Br-\-\) 


(ho{l3  +  2Co)  -  Qi(cos  Bi  -  1))  +  O 


(16) 
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It  turns  out,  however,  that  for  larger  amplitude  waves  this  expansion  is  one  of 
the  most  significant  error  sources.  It  is  therefore  useful  to  carry  the  expansion  of 
(15)  to  second  order  which  yields  a  quadratic  equation  in  Again  we  can  solve 

for  and  eliminating  the  false  root  we  then  get  the  second  order  expression  for 

Qx,r- 

^  cos  9r  (coho{cos  +  1)  +  Qi(cOS  6i  -  cos  Or  -  §)) 

^  (2cos9,  +  |) 


(4cos^r  +  3)((^f(|  -  cosOj)  +  QiCoho{cosOi  -  1)  -  CphKP  +  2Co))  ^  ^ 
^  {CoK{cOs6r  +  1)  +  Qi{cos0i  -  cos  Or  -  1))^ 


(17) 


These  equations  have  two  unknowns,  Qx.r  and  Or,  which  can  be  determined  by 
realizing  that 

Or  =  arctan 

\^x,r  J 

The  additional  unknown  can  be  solved  by  using  (9)  which  is  rewritten  as; 

Qy,r  =  7(^0  +  0  ~  Qy>i 

in  which  Qy^i  is  specified  and  q  is  determined  by  integration  of  the  last  of  the 
characteristic  equations  (7). 

In  these  expressions  is  the  Riemann-variable  updated  to  the  next  time  level 
by  (5).  Qi  is  the  total  flux  of  the  known  incoming  wave  at  the  same  time  level.  From 
(16)  or  (17),  and  (18)  and  (19)  we  can  find  the  unknowns  Q,r,r  and  Or  iteratively. 
With  the  incoming  wave  known  through  specification,  the  boundary  value  of  total 
flux  Qj;  can  determined  at  the  next  time  step. 


4  Reflection  properties 


The  absorption  properties  of  the  boundary  condition  are  tested  for  a  unidirec¬ 
tional  wave  in  a  domain  of  constant  depth  for  various  angles  of  incidence  and  wave 
amplitudes.  In  the  following  example  waves  are  generated  at  the  x  =  0  and  y  =  0 
boundaries  and  absorbed  at  all  four  boundaries  for  both  versions  of  the  bound¬ 
ary  condition.  The  physical  parameters  are  still  water  depth  ho  and  wave  length 
A  =  100  ho,  which  yields  a  period  T  =  lOO^h^.  The  numerical  parameters  are 
Ax  =  Ay  =  A/60  and  =  T/lOO  which  yields  a  Courant  number  of  0.6.  A 
second-order  predictor-corrector  numerical  scheme  is  used  for  all  tests. 

Linear  waves 

In  the  first  test,  sinusoidal  waves  with  a  small  amplitude  A/ ho  =  0.01  are  prop¬ 
agated  using  the  linear  equations  and  absorbed  using  the  boundary  condition  that 
applies  the  lowest  order  expansion  of  Qx,r,  (16).  Sinusoidal  waves  are  specified  as 
the  initial  condition.  This  case  was  previously  shown  in  Van  Dongeren  et  al.  (1994). 

The  reflection  properties  are  computed  for  a  square  domain 

ni-{(a;,j/);0<x<  A,0<y<A}  (20) 

where  we  want  to  make  sure  that  the  reflections  are  caused  by  one  absorbing  bound¬ 
ary  (at  a;  =  A)  only.  This  is  accomplished  as  follows.  Solutions  are  computed  in  two 
domains;  a  rectangular  domain 

02  =  {(x,y):0<a;<A,0<?/<3A}  (21) 


and  in  a  larger,  square  domain. 


see  Fig.  4  for  a  definition  sketch. 


The  non-generating  boundaries  in  domain  fla  are  placed  so  far  away  that  they 
have  no  effect  on  the  solution  in  the  smaller  domain  Oi  during  the  duration  of 
the  simulation.  Therefore  in  the  smaller  domain  we  can  consider  the  fis-solution 
free  of  reflection  errors.  Similarly,  in  domain  ^2  the  non-generating  boundary  at 
y  =  3  A  will  not  influence  the  solution  in  the  smaller  domain  .  Hence  the  difference 
between  the  two  solutions  can  only  be  caused  by  the  absorbing  boundary  at  x  =  A. 
The  two  solutions  are  subtracted  from  each  other  at  the  instant  in  time  =  T j  cos  9i 
when  the  initial  condition  has  propagated  out  of  Hi  and  the  difference  is  normalized 
by  the  amplitude  A: 


Ci{x,  y,t^-,9i) 


A 


(23) 


where  and  are  the  solutions  for  the  test  runs  in  the  domains  0,2  and  H3, 
respectively. 


Eq.  (23)  yields  a  spatial  picture  of  the  reflection  error  in  Hi  due  to  the  absorbing 
boundary  condition  at  x  =  A.  Figs.  5a,  c  and  e  show  contours  of  the  spatial  errors 
for  three  angles  of  incidence:  0,-  =  0,  |  and  f  where  9i  is  defined  as  the  angle  between 
the  direction  of  propagation  and  the  x-axis.  The  errors  are  of  the  order  0.005  or 
0.5%. 


Note  that  the  boundary  condition  is  derived  based  on  the  nonlinear  equations  in 
characteristic  form  while  the  waves  themselves  are  run  under  the  linear  equations, 
which  is  in  itself  inconsistent  but  allowable  for  small  amplitude  waves. 


Nonlinear  waves 


To  show  the  effect  of  the  nonlinear  terms  in  the  equations,  a  second  test  is  run  for 


the  same  parameters  and  for  the  same  angles  of  incidence  but  using  the  nonlinear 
equations.  The  waves  generated  at  the  boundary  are  again  sinusoidal.  The  spatial 
errors  are  plotted  in  Figs.  5b,  d  and  f.  Comparison  to  the  previous  case  shows  that 
including  the  nonlinear  terms  in  the  governing  equations  increases  the  error  by  a 
factor  2:  they  are  now  of  the  order  of  1%.  This  is  due  to  the  fact  that  in  the  first 
order  approximation  (16),  Eq.  (14)  reduces  to  the  linear  relationships  of  constant 
form 

Qx,i  ~  Ci  COS  Qx,t  ~  Cq  COS  9^  (^4) 

where  is  given  by  (12).  This  means  that  the  first  order  boundary  condition  absorbs 
constant  form  (linear)  waves  better  than  waves  that  change  shape. 

As  can  be  seen  in  Fig.  5,  the  errors  are  spatially  dependent.  In  order  to  obtain 
a  single  measure  of  the  error  as  a  function  of  the  angle  of  incidence,  the  T^-error  is 
computed  (Strikwerda,  1989).  The  T^-error  is  defined  as  the  squared  difference  of 
(qz  and  Coa  evaluated  in  the  domain  fli  and  normalized  by  the  RMS  of  the  larger 
area  solution  at  the  time  instant  t  =  when  the  initial  condition  has  propagated 
out  (which  is  different  for  each  9i): 

, .  „ ,  \/En, (&, 

=  - - ,  - 

Figure  6  shows  L^-error  for  the  range  of  angles  of  incidence  of  6i  =  [0,  |,  ...  ,  |] 

for  both  the  linear  and  nonlinear  low-amplitude  cases  described  above.  Also  plotted 
in  Figure  6  is  the  error  incurred  when  the  Sommerfeld  radiation  condition  (2)  is 
applied  for  sinusoidal  waves.  This  condition  shows  near-perfect  absorption  for  waves 
of  normal  incidence  but  shows  large  errors  for  more  obliquely  incident  waves.  In  fact, 
the  error  is  100%  for  glancing  angles.  In  contrast,  the  errors  due  to  the  boundary 


condition  derived  in  Section  3  are  of  the  order  of  0.5%  to  1%  for  the  whole  range  of 
angles  of  incidence,  which  is  acceptable  for  most  applications. 

However,  the  error  is  a  function  of  the  nonlinearity  parameter,  8  =  Ajho-  It  can 
be  shown  that  if  (16)  is  applied  in  the  boundary  condition,  the  error  e2  ~  5.  In  a 
third  test,  the  model  is  run  for  the  same  parameters  as  the  previous  test  but  with 
a  wave  with  a  ten  times  larger  amplitude,  Ajho  =  0.1.  The  spatial  variation  of 
the  reflection  errors,  see  Fig.  7a,  c  and  e,  are  about  one  order  of  magnitude  larger 
than  in  the  previous  test,  as  one  should  expect.  Fig.  8  (dashed  line)  shows  that  the 
L^-error  versus  9i  for  a  single  absorption  boundary  is  about  10%,  or  6  *  100%,  which 
is  too  large  for  practical  purposes. 

As  mentioned  in  Section  3,  a  major  error  source  for  large  amplitude  waves  is  the 
first  order  expansion  of  (15).  For  larger  amplitudes,  it  is  therefore  advantageous  to 
use  the  second  order  approximation  (17)  as  the  boundary  condition.  The  spatial 
errors  as  calculated  by  (23)  are  shown  in  Fig.  7b,  d  and  f.  We  see  that  the  error 
is  now  reduced  by  a  factor  5  compared  to  the  left-hand  side  panels.  For  a  single 
absorbing  boundary,  the  T^-error  versus  the  angle  of  incidence  di  as  calculated  by 
(25)  is  shown  in  Fig.  8  (solid  line).  It  has  a  magnitude  of  about  3%,  which  is  close  to 
the  theoretical  error  e2  ~  6/4  for  the  second  order  approximation.  For  our  purposes, 
this  magnitude  of  the  error  for  medium  height  waves  is  acceptable. 

Since  the  second  order  approximation  (17)  is  explicit  in  Qx,ri  it  does  not  require 
considerably  more  computational  time  than  the  linear  expression  (16)  and  it  will 
therefore  be  used  in  the  remainder  of  this  paper. 
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5  Example  of  simultaneous  absorption  and  gen¬ 
eration 


To  illustrate  the  application  of  the  proposed  boundary  condition  for  the  case 
of  simultaneous  absorption  and  generation  of  waves  at  one  boundary,  consider  a 
domain  with  an  absorbing-generating  boundary  at  the  a:  =  0  boundary  and  a  wall 
at  a;  =  4.25A.  From  a  cold  start,  incoming  waves  are  generated  at  normal  incidence 
from  t  =  0  T  -  19  r,  tapered  with  a  hyperbolic  tangent  function  during  the  first  and 
last  period  of  generation  in  order  to  eliminate  transients  due  to  shocks.  This  wave 
train  will  reflect  off  the  wall  at  a;  =  4.25A  and  produce  a  standing  wave  until  the 
incident  waves  are  turned  off  at  T  =  19T.  Parameters  used  are:  water  depth  ho, 
wave  length  A  =  lOOho,  T  =  X  j  yjg  ho,  A  =  O.Ol/io,  Aa:  =  Ay  =  ^  and  Cr  =  0.6. 
The  numerical  scheme  for  the  linearized  equations  is  used  in  this  example. 

Fig.  9a  shows  the  time  series  of  the  specified  volume  flux  of  the  incident  waves  at 
a;  =  0,  normalized  by  Co  ho-  The  effect  of  the  hyperbolic  tangent  function  is  clearly 
visible  as  the  amplitude  grows  to  its  full  value  in  little  more  than  one  wave  period. 
Fig.  9b  shows  the  flux  of  the  outgoing  wave  train  at  the  same  point  as  calculated  by 
(17).  The  outgoing  wave  train  is  a  near-perfect  mirror  image  of  the  incident  wave 
train  except  for  some  small  trailing  waves  around  t  =  28T.  Still  water  is  recovered 
almost  immediately  after  the  outgoing  wave  passes  through  the  a;  =  0  boundary, 
which  indicates  that  the  reflections  are  small.  The  time  series  of  the  total  flux  in 
Fig.  9c,  which  is  the  sum  of  the  two  time  series  above,  shows  first  a  progressive 
wave  in  the  -\-x  direction,  then  the  anti-node  of  the  flux  of  a  standing  wave,  then  a 
progressive  wave  in  the  -x  direction  and  finally  still  water.  This  Figure  qualitatively 
shows  that  the  absorbing-generating  boundary  condition  works  very  well. 


In  order  to  measure  the  reflection  error  due  to  the  absorbing-generating  boundary 
condition,  the  envelope  of  the  standing  wave,  (,  is  calculated- when  the  front  of  the 
incident  wave  train  has  propagated  through  the  domain  and  has  almost  reached  the 
a:  =  0  boundary  (which  happens  at  t  =  8.5  T).  The  envelope  is  shown  in  Fig.  10a 
(solid  line).  This  Figure  does  not  show  the  entire  domain  because  the  tapered  front 
of  the  wave  does  not  produce  a  standing  wave  at  this  time.  The  standing  wave  is 
a  result  of  the  summation  of  the  incident  wave  train  and  its  reflection  off  the  wall 
and  has  a  maximum  amplitude  of  2  A.  When  the  front  of  the  wave  train  reaches  the 
open  boundary  at  t  —  8.5  T,  it  will  be  absorbed  almost  perfectly  as  indicated  in  Fig. 
9b.  However,  a  small  portion  of  the  wave  with  amplitude  will  be  reflected  and 
propagate  back  in  the  +x  direction.  This  small  error  wave  will  itself  reflect  from  the 
back  wall  and  produce  a  standing  wave  of  its  own.  Again,  the  envelope  of  the  total 
standing  wave  (the  standing  wave  due  to  the  specified  wave  and  the  error  wave)  can 
be  calculated  at  a  time  just  before  this  small  error  wave  reaches  the  open  boundary, 
which  happens  at  t  =  17  T.  This  envelope,  C  +  Co  has  a  maximum  amplitude  of 
2{A  +  Ae)  and  is  shown  as  the  dashed  line  in  Fig.  10a.  Due  to  the  smallness  of  the 
error  wave,  the  dashed  line  is  almost  indistinguisable  from  the  solid  line.  A  measure 
of  the  error  in  amplitude  caused  by  the  generating- absorbing  boundary  can  then  be 
defined  as  the  difference  between  the  two  envelopes  normalized  by  the  amplitude  of 
the  original  standing  wave  or 


(26) 


which  is  shown  in  Fig.  10b  and  has  a  maximum  of  about  0.2%. 
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6  Discussion 


The  numerical  tests  described  above  show  that  the  reflection  errors  due  to  the 
present  boundary  condition  are  of  the  order  of  only  a  few  percent  for  cases  of 
absorption-generation  as  well  as  absorption  only.  This  is  a  remarkable  improve¬ 
ment  over  the  absorption  properties  of  the  widely-used  radiation  conditions  based 
on  the  wave  equation,  which  only  absorb  waves  at  one  specific  angle  of  incidence 
perfectly  and  show  large  errors  for  other  angles  of  incidence. 

In  its  present  form  the  boundary  condition  has  one  major  drawback.  It  cannot 
absorb  multiple  waves  with  a  difference  between  angles  of  incidence  larger  than  |. 
This  is  due  to  the  fact  that  by  using  (14)  the  set  of  multiple  waves  at  the  boundary 
is  essentially  approximated  by  one  representative  progressive  wave  which  is  not  valid 
when  two  or  more  wave  trains  intersect  at  oblique  angles. 

It  is  emphasized  that  for  simplicity  the  present  form  of  the  boundary  condition 
does  not  account  for  time-varying  or  steady  currents.  However,  the  mathematical 
modification  of  either  (16)  or  (17)  to  include  currents  is  straightforward.  The  real 
problem,  as  stated  in  Section  2,  lies  in  the  philosophical  distinction  between  currents 
and  long  waves  and  in  the  fact  that  currents  would  have  to  be  known  a  priori  and 
specified  along  all  boundaries  of  the  domain. 

7  Conclusions 

In  this  paper  two  orders  of  an  absorbing-generation  boundary  condition  for  the 
nonlinear  shallow  water  equations  are  derived  based  on  the  Method  of  Characteris¬ 
tics.  Numerical  tests  show  that  by  using  this  boundary  condition,  reflection  errors 


can  be  limited  to  a  few  percent  of  the  incident  wave  amplitude  for  the  full  range  of 
angles  of  incidence,  which  is  an  improvement  over  the  classical  radiation  conditions. 
Unlike  those  radiation  conditions,  the  present  boundary  condition  allows  simulta¬ 
neous  specification  of  an  incident  wave  train  and  absorption  of  an  outgoing  wave 
train  at  the  same  boundary,  which  makes  it  particularly  suitable  for  application  on 
artificial  oceanside  boundaries  in  shallow  water  models. 


Acknowledgments 


This  work  is  a  result  of  research  sponsored  by  NOAA  Office  of  Sea  Grant,  De¬ 
partment  of  Commerce,  under  Award  No.  NA  56  RG  0147  (Project  No.  R/OE-17) 
and  by  the  U.S.  Army  Research  Office,  University  Research  Initiative  under  Con¬ 
tract  No.  DAAL03-92-G-0016.  The  U.S.  Government  is  authorized  to  produce  and 
distribute  reprints  for  government  purposes  notwithstanding  any  copyright  notation 
that  may  appear  herein. 


Appendix  I.  References 

Abbott,  M.B.  (1979).  Computational  Hydraulics.  Elements  of  the  Theory  of  Free 
Surface  Flow.  Pitman  Publishing,  324  p. 

Broeze,  J.  and  E.F.G.  Van  Daalen  (1992).  Radiation  boundary  conditions  for  the 
two-dimensional  wave  equation  from  a  variational  principle.  Math.  Comp., 
vol.  58,  no.  197,  pp.  73-82. 

Engquist,  B.  and  A.  Majda  (1977).  Absorbing  boundary  conditions  for  the  numer¬ 
ical  simulation  of  waves.  Math.  Comp.,  vol.  31,  no.  139,  pp.  629-651. 


19 


Givoli,  D.  (1991).  Non- reflective  boundary  conditions.  J.  of  Comp.  Physics.,  vol. 
94,  pp.  1-29. 

Hibberd,  S.  (1977).  Surf  and  run-up.  Ph.D.  dissertation,  School  of  Mathematics, 
University  of  Bristol,  124  p. 

Higdon,  R.L.  (1986).  Absorbing  boundary  conditions  for  difference  approximations 
to  the  multi-dimensional  wave  equation.  Math.  Comp.,  vol.  ^7,  no.  176,  pp. 
437-459. 

Higdon,  R.L.  (1987).  Numerical  absorbing  boundary  conditions  for  the  wave  equa¬ 
tion.  Math.  Comp.,  vol.  46,  no.  179,  pp.  65-90. 

Keys,  R.G.  (1985).  Absorbing  boundary  conditions  for  acoustic  media.  Geophysics, 
vol.  50,  no.  6,  pp.  892-902. 

Kobayashi,  N.,  A.K.  Otta  and  I.  Roy  (1987).  Wave  reflection  and  run-up  on  rough 
slopes.  Journal  of  Waterway,  Port,  Coastal  and  Ocean  Engineering,  ASCE, 
vol.  ns,  no.  3,  pp.  282-298. 

Sommerfeld,  A.  (1964).  Lectures  on  Theoretical  Physics.  Academic  Press,  New 
York. 

Strikwerda,  J.C.  (1989).  Finite  difference  schemes  and  partial  differential  equa¬ 
tions.  Wadsworth  k  Brooks/Cole  Mathematical  Series.  Wadsworth,  Belmont, 
California,  386  pp. 

Svendsen,  LA.  (1974).  Cnoidal  waves  over  a  gently  sloping  bottom.  5enes  pa¬ 
per  #  6,  Institute  of  Hydrodynamics  and  Hydraulic  Engineering,  Technical 
University  of  Denmark,  181  pp. 


20 


Svendsen,  I.A.  and  P.  Justesen  (1984).  Forces  on  slender  cylinders  from  very  high 
waves  and  spilling  breakers.  Symposium  on  description  and  modelling  of  di¬ 
rectional  seas.  Copenhagen  1984.  Chapter  D-7,  16  pp. 

Svendsen,  I.A.  and  U.  Putrevu  (1995).  Surf-zone  hydrodynamics.  Research  report 
CACR-95-02,  University  of  Delaware,  71  pp. 

Van  Dongeren,  A.R.,  F.E.  Sancho,  I.A.  Svendsen  and  U.  Putrevu  (1994).  SHORE- 
CIRC:  A  quasi  3-D  nearshore  model.  Proc.  of  the  24th  ICCE,  pp.  2741-2754. 

Van  Dongeren,  A.R.,  I.A.  Svendsen  and  F.E.  Sancho  (1995).  Application  of  the 
Q-3D  SHORECIRC  model  to  surfbeat.  Proc.  Coastal  Dynamics.,  pp.  233-244. 

Verboom,  G.K.,  G.S.  Stelling  and  M.J.  Oflicier  (1981).  Boundary  conditions  for  the 
shallow  water  equations.  In:  Abbott,  M.B.  and  J.A.  Cung  Eds..,  Engineering 
Applications  of  Computational  Hydraulics,  vol.  I,  pp.  230-262. 

Verboom,  G.K  and  A.  Slob  (1984).  Weakly-reflective  boundary  conditions  for  two- 
dimensional  shallow  water  flow  problems.  Adv.  Water  Resources,  vol.  7,  pp. 
192-197. 

Whitham,  G.B.  (1974).  Linear  and  nonlinear  waves.  John  Wiley  and  Sons,  New 
York,  636  pp. 


21 


Appendix  II.  Notation 


The  following  symhoi 

[5  are  used  in  this  paper: 

A 

- 

wave  amplitude 

Ae 

- 

amplitude  of  reflected  wave 

c 

- 

nonlinear  shallow  water  phase  speed 

Co 

- 

linear  shallow  water  phase  speed 

f.Jy 

- 

forcing  terms  in  NSW-equations 

Fp+,Fp-,F^ 

- 

forcing  terms  in  Equations  in  Characteristic  Form 

9 

- 

gravity 

h 

total  water  depth 

ho 

- 

still  water  depth 

P 

- 

eigenvector  space 

Q 

- 

total  flux 

Qx^  Qy 

- 

total  flux  in  x  and  y  directions 

Qi 

- 

total  flux  of  the  incoming  wave 

Qx,ii  Qy,i 

- 

flux  of  the  incoming  wave  in  x  and  y  directions 

QxjT^  Qy,T 

- 

flux  of  the  outgoing  wave  in  x  and  y  directions 

Q 

- 

total  flux  of  the  current 

U,  V 

- 

shortwave  averaged  and  depth  averaged  horizontal  velocities 

T 

- 

wave  period 

- 

Riemann  variables 

6 

- 

nonlinearity  parameter 

Cl,  C2,  £3 

- 

reflection  errors 

_c_ 

- 

shortwave  averaged  surface  elevation 

_Cii  Cj 

- 

surface  elevation  of  incoming  and  outgoing  wave 

C02  ,  Cfis 

- 

surface  elevation  computed  in  the  domains  ^2  and  O3 

c 

- 

infragravity  wave-averaged  surface  elevation 

c,c 

- 

envelopes  of  the  standing  wave 

A 

- 

wave  length 

On 

- 

angle  of  incidence  of  the  radiated  waves 

Oi,  Or 

- 

angle  of  incidence  of  incoming  and  outgoing  wave 

d 

- 

angle  between  boundary  and  a;-coordinate  axis 

)  ^2 ,  ^3 

- 

domains  used  in  reflection  tests 
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Figure  3:  Definition  sketch  of  the  characteristics. 
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Figure  5:  Reflection  errors  vs  (a:/A,  y/A)  for  wave  amplitude  A//io  =  0.01  using  (16) 
(a)  linear  scheme,  di  =  0,  at  =  T;  (6)  nonlinear  scheme,  Bi  =  0,  at  =  T; 
(c)  linear  scheme,  Bi  =  |,  at  T  =  T/cos|;  (d)  nonlinear  scheme,  Bi  -  |,  at 
r  =  r/cos|;  (e)  linear  scheme,  Bi  =  |,  at  T  =  T/ cos  f ;  (/)  nonlinear  scheme, 
0i  =  |,  att’^  =  r/cosf. 
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Figure  7:  Reflection  errors  vs  (x/A,j//A)  for  wave  amplitude  kjho  =  0.1  and  a 
nonlinear  scheme  (a)  Oi  =  0,  at  =  T,  first  order  BC;  (6)  Oi  =  0,  at  =  T, 
second  order  BC;  (c)  Oi  =  at  T  =  r/cos|,  first  order  BC;  (d)  Oi  =  at 
r  =  r/cos|,  second  order  BC;  (e)  Oi  =  at  =  r/cos|,  first  order  BC;  (i) 
Oi  —  aX  —  T/ cos  f,  second  order  BC. 


Figure  9:  Time  series  at  a;  =  0:  (a)  incident  wave  flux,  (b)  outgoing  wave  flux,  (c) 
total  flux. 


